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PREFACE 


This  final  technical  report  was  prepared  by  the  Texas  A&M  Research 
Foundation  under  Grants  DAHC04-74-G-0184  (1  June  1974  - 10  June  1976) 
and  DAAG29-76-G-0241  (11  June  1976  - 31  December  1977)  . These  grants 
were  the  follow  on  efforts  after  Project  Themis. 

In  the  summer  of  1968,  the  U.S.  Department  of  Defense  through 
the  U.S.  Ariny  Research  Office  at  Durham  sponsored  a program  of  research  on 
Aircraft  Dynamics  for  Subsonic  Flight  at  Texas  A&M  University.  The  work 
was  funded  under  Project  Themis  (Contract  DAHC04-69-C-0015)  which  had 
been  established  a year  earlier  in  response  to  the  late  President  Johnson's 
request  that  each  Federal  Agency  should  help  to  develop  new  centers  of 
excellence  in  areas  relevant  to  its  goals.  By  means  of  Project  Themis, 
the  Department  of  Defense  hoped  to  meet  in  part  its  long-term  research 
needs,  strengthen  more  of  the  nation's  universities,  increase  the  number 
of  institutions  performing  research  of  high  quality  and  achieve  a wider 
distribution  of  research  funds.  In  this. way,  it  planned  to  enhance  the 
United  States'  academic  capabilities  in  science  and  technology. 

After  Project  Themis  expired  in  May  1973,  Army  Research  Office 
awarded  another  Grant  DA-ARO-D-31-124-71  G153  which  was  followed  by  the 
present  grants.  Dr.  James  J.  Murray  was  the  Technical  Monitor  for  Project 
Themis  efforts.  Dr.  Robert  E.  Singleton  was  the  Technical  Monitor  for  the 
present  effort.  Dr.  W.  P.  Jones  was  the  Program  Manager  for  Project  Themis, 
and  later,  was  one  of  the  Co-Principal  Investigators  along  with  the  senior 
author  of  this  report  till  he  retired  from  Texas  A&M  University  in  August 
1975.  At  that  time,  B.  M.  Rao  became  the  Principal  Investigator  for  the 
grants. 

Texas  A&M  University  benefited  greatly  from  this  entire  program  sponsored 
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by  the  Army  Research  Office  to  support  our  research  work  with  a view  to 
developing  a center  of  excellence  for  the  study  of  aerodynamic  and  dynamic 
problems  of  aircrafts  and  helicopters.  Before  this  program,  the  University 
did  not  have  an  active  research  participation  in  this  field.  During  this 
period,  more  than  twenty  students  who  worked  in  these  research  projects 
received  their  graduate  degrees  (Master's  and  Ph.D.'s).  Several  of  these 
students  are  presently  working  for  the  Aerospace  industry  and  the  U.  S. 
Government.  During  the  same  period,  some  of  these  students  won  awards  for 
their  outstanding  research  papers  in  AIAA  Annual  Regional  and  National 
Paper  Competitions.  This  long-term  funding  from  the  Army  Research  Office  has 
had  a two-fold  effect:  first,  some  basic  aerodynamic,  efficient  compu- 
tational programs  have  been  developed,  and  are  being  used  by  scientists 
in  the  field.  Secondly,  several  students  have  received  excellent  edu- 
cational training  and  research  experience  working  in  the  sponsored  research 
projects  under  professional  engineers. 

The  authors  gratefully  acknowledge  Dr.  Robert  E.  Singleton  for  his 
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sustained  interest  in  the  program,  and  for  his  encouragement  and  cooperation 
throughout  the  present  grant  period. 
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I.  INTRODUCTION 


Although  man  has  been  flying  aircrafts  for  many  years,  a sizeable 
range  of  problems  continue  to  exist.  This  is  due  partly  to  the  appearance 
of  new  vehicles  such  as  V/STOL  aircraft  and  partly  to  the  changes  in  the 
aircrafts'  operational  requirements.  Hazen^  emphasized  the  need  for  de- 
voting greater  effort  to  the  study  of  subsonic  aerodynamics.  He  noted 
that  most  current  problems  are  dynamic  in  nature,  and  are  concerned  with 
aircraft  maneuverability,  response  to  atmospheric  turbulence,  aeroelastic 
phenomena  and  non-linear  aerodynamic  problems  arising  from  separated 
flows.  A better  understanding  of  these  problems  is  hinged  upon  the 
availability  of  efficient  and  accurate  airloads  prediction  techniques  for 
lifting  surfaces.  Comprehensive  literature  survey  for  the  steady  subsonic 
lifting  line  and  lifting  surface  theories  is  presented  in  reference  2,  and 
a brief  summary  of  this  survey  follows. 

One  most  basic  approach  in  predicting  spanwise  loads  on  subsonic, 
three-dimensional  wings  is  the  lifting  line  technique  due  to  Prandtl^ 
which  is  applicable  to  large  aspect  ratio  unswept  planforms.  The  wing  is 
represented  by  a bound  vortex  at  the  quarter  chord,  and  a trailing  vortex 
sheet  extending  from  the  bound  vortex  to  infinity.  The  analysis  is  based 
on  the  assumption  that  at  every  point  along  the  span  the  flow  is  essentially 
two-dimensional:  the  induced  angle  of  attack  produced  by  the  trailing 
vortex  wake  is  used  to  correct  the  two-dimensional  lift  of  the  bound  vortex 
to  account  for  the  finite  wing.  The  Prandtl  lifting  line  concept  does  not 
produce  the  chordwise  pressure  distribution  and  is  certainly  limited  by 
restrictions  on  wing  planform. 

Weissinger^  extended  Prandtl 's  lifting  line  concept  to  swept  wings  by 
modifying  the  straight  bound  vortex  to  a V-shaped  vortex  at  the  quarter 
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chord  of  the  wing.  In  this  technique,  known  as  the  Weissinger-L  method, 
the  velocity  induced  at  the  three-quarter  chord  by  the  trailing  wake 
system  and  the  inclined  bound  vortex  pair  is  equated  to  the  component  of 
free  stream  velocity  normal  to  the  wing  surface.  This  specification  is 
justified  by  the  fact  that  in  two-dimensional  theory  this  arrangement 
yields  the  correct  circulation.  An  interesting  variation  of  the  Weissinger- 
L approach  is  a finite  element  method  by  Campbell'*,  in  which  the  swept 
bound  vortex  pair  and  the  trailing  wake  are  approximated  by  a system  of 
rectangular  horseshoe  vortices  arranged  with  their  midpoints  along  the 
quarter  chord.  This  approach  is  extended  by  Blackwell^  to  wing-tail 
combinations  and  to  wings  with  pylons  and  end  plates. 

Although  the  lifting  line  methods  predict  the  spanwise  loading  ac- 
curately for  wings  of  conventional  planforms,  the  extremes  of  sweep,  taper 
and  aspect  ratios  encountered  in  modern  lifting  surface  design,  along  with 
a need  for  a complete  description  of  the  surface  load  distribution  have 
resulted  in  the  development  of  the  lifting  surface  theories.  The  lifting 
surface  problem  has  been  approached  approximately  by  multiple  lifting-line 
and  discrete  loading  element  schemes.  The  direct  approach^ employs 
a series  substitution  in  the  governing  integral  equations  relating  the  normal 
velocities  on  the  surface  to  the  pressure  load,  or  an  equivalent  quantity 
across  the  surface.  The  series  is  derived  from  subsonic  thin  airfoil 
theory  for  chordwise  loading  and  from  the  lifting  line  approach  for  spanwise 
loading.  The  series  substitution  allows  evaluation  of  the  lifting  surface 
integral  equation  at  a number  of  surface  points  equal  to  the  desired  number 
of  series  coefficients,  or  for  a larger  number  of  surface  points  in  u 
'least  squares'  sense. 

Falkner^  conceived  the  discrete  loading  element  scheme  in  which  the 
lifting  surface  is  replaced  by  a system  of  concentrated  horseshoe  vortices. 


and  their  strengths  computed  by  satisfying  the  tangential  flow  condition 
at  surface  points  equal  in  number  to  the  unknown  circulations.  High 
speed  digital  computers  have  enabled  investigators  to  refine  the  discrete 
loading  element  concept  to  allow  solutions  for  larger  systems  of  complex 
geometries.  Woodward^®  represents  an  arbitrary  wing-body  combination  in 
steady  motion  by  a system  of  source  doublet  and  vortex  singularities. 

Lift  and  volume  effects  of  the  body  are  effected  by  line  singularities 
along  the  body  axis,  while  wing  lift  and  wing-body  interference  are  rep- 
resented by  planar  singularities  over  a finite  number  of  quadrilateral 
panels  on  the  wing  and  interface  region  of  the  body.  The  strengths  of 
the  singularities  are  approximated  as  constant  over  each  panel  and  de- 
termined so  as  to  satisfy  flow  tangency. 

The  problems  of  unsteady  flows  have  been  studied  by  several  in- 
vestigators. Analytical /numerical  solutions  of  the  governing  equations 
of  motion  are  obtained  by  requiring  that  these  solutions  satisfy  the 
appropriate  boundary  conditions.  A brief  literature  survey  is  presented 

• t 

in  reference  11.  Recently,  discrete  or  finite  element  schemes  have 
been  developed  for  oscillatory  subsonic  flows  with  good  success,  but 
these  schemes  include  several  approximations  to  facilitate  solution  of  the 
complex  mathematical  formulations.  Albano  and  Rodden^^,  and  Kalman 
et  al^^  f,ave  demonstrated  a wide  variety  of  applications  of  the  doublet- 
lattice  method  which  is  an  extension  of  the  vortex- lattice  method  of 
Falkner  for  unsteady  flows.  In  this  technique,  the  lifting  surface  to 
be  analyzed  is  divided  into  a large  number  of  planar  panels.  Each  panel 
is  then  replaced  by  a horseshoe  vortex  representing  the  steady  portion 
of  the  flow,  and  an  oscillatory  doublet  distribution  representing  the 
unsteady  part  of  the  flow.  The  collocation  method  is  then  used  to  solve 
for  the  vortex  and  the  doublet  strengths  incorporating  the  appropriate 
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boundary  conditions.  Reference  12  includes  the  analysis  of  an  oscillating 
airfoil,  an  airfoil  in  an  oscillating  flow,  and  variable  geometry  wings 
with  and  without  flaps  while  in  reference  13,  the  doublet  lattice  method 
has  been  used  to  analyze  non-planar  configurations  such  as  T-tail,  a wing 
with  end  plates,  and  annular  wings.  This  method  is  very  restricted  in  its 
applications  to  non-planar  configurations  such  as  wing-body  combinations. 

As  each  doublet  panel  has  associated  with  it  a trailing  wake  which  extends 
infinitely  aftward  in  a direction  parallel  to  the  free  stream,  the 
placement  of  panels  on  a body  cannot  be  such  as  would  allow  the  trailing 
wake  of  any  panel  to  penetrate  the  body,  for,  this  would  constitute  a 
physically  unrealistic  situation.  Thus  a body  representation  is  generally 
restricted  to  being  an  annular  wing  of  constant  cross  section  without  any 
inclination  to  the  flow. 

The  analysis  of  a cascade  of  blades  in  unsteady  flows  has  played  an 
important  role  in  compressor  and  turbine  blade  design.  Chang  and  Chu^‘* 
have  analyzed  a cascade  of  blades  in  synchronized  oscillations.  Their 
solution  involved  a modified  form  of  The'odoi'sen's  function  which  reduced 
to  Theodorsen's  function  as  the  gap  between  two  adjacent  airfoi Is  became 
infinite.  Another  approach  by  Schorr  and  Reddy^^  has  been  used  to  analyze 
the  oscillatory  flow  over  cascades  with  varying  solidity  and  stagger  angle. 

In  their  analysis,  the  airfoils  in  the  cascade  are  replaced  by  a distribution 

i 

of  vortices  and  source-sink  pairs  whose  strengths  are  determined  for  the  j 

appropriate  boundary  conditions.  Since  the  formulation  is  very  complicated, 
they  are  able  to  obtain  only  an  approximate  solution  by  replacing  the 
kernel  function  with  a polynomial. 

Jones  and  Moore^^,  and  Rao  and  Jones^^  solved  the  problems  of  a 
cascade  of  oscillating  airfoils  in  incompressible  flow  and  oscillating 
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incompressible  flow  over  a cascade  of  airfoils,  respectively.  They 
used  the  present  lifting  surface  method  based  upon  the  velocity  potential 
formulation,  and  were  able  to  develop  a very  efficient  computational 
scheme  as  an  alternative  to  the  existing  lifting  surface  methods  based  on 
acceleration  potential  formulation. 

The  problem  of  wind-tunnel  wall  interference  effects  on  oscillating 
airfoils  in  subsonic  flow  constitutes  a special  case  of  the  blade-row 
problem,  namely,  the  case  in  which  stagger  angle  is  zero  and  adjacent 
blade  angles  are  180°  out  of  phase.  Runyan  and  Watkins treated  the 
effect  of  wind-tunnel  walls  on  the  aerodynamic  forces  of  an  oscillating 
airfoil  in  a subsonic  compressible  flow.  They  simulated  the  walls  by 
placing  images  at  appropriate  distances  above  and  below  the  wings  such 
that  the  condition  of  zero  normal  velocity  at  the  tunnel  walls  is  satisfied. 
In  their  formulation,  the  integral  equation  between  the  downwash  and  the 
oscillatory  doublet  strength  involved  a series  of  Hankel  functions.  The 
solution  of  this  resulting  equation  yields  the  oscillatory  doublet  strength 

' I 

as  a function  of  chordwise  location  on  the  airfoil  and  these  doublet 
strengths  were  used  to  obtain  acceleration  potentials  and  hence  the  pressure 
distribution  on  the  airfoil.  An  important  result  shown  was  that,  for  certain 
conditions  of  frequency,  tunnel  height,  Mach  number,  the  tunnel  and  wing 
may  form  a resonant  system  so  that  the  forces  on  the  wing  are  greatly 
changed  from  the  condition  of  no  tunnel  walls. 

Jones*®,  independently  developed  a theory  for  estimating  the  effect 
of  wind-tunnel  walls  on  the  airforces  acting  on  an  airfoil  oscillating  in  a 
subsonic  airstream.  In  his  formulation,  the  integral  equation  between  the 
downwash  and  the  doublet  strength  (modified  velocity  potential)  was  es- 
sentially similar  to  the  equation  developed  by  Runyan  and  Watkins.  Jones 
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suitined  the  series  of  Hankel  functions  in  the  integral  and  developed  a numerical 
procedure  applicable  to  the  low  frequency  values.  The  results  obtained  for 
a special  case  compared  well  with  available  experimental  results. 

Lane  and  Friedman^o  developed  a method  for  calculating  the  aerodynamic 
forces  experienced  by  a cascade  with  arbitrary  stagger  and  interblade 
phase-lag  angles.  Their  method  utilized  Fourier  transforms  of  blade 
pressure- jump  functions  which  permit  the  kernel  function,  appearing  in  the 
integral  equation  relating  the  upwash  to  the  pressure-jump  transform,  to  be 
expressed  in  a closed  form,  by  assuming  a commonly  used  trigonometric  series 
expansion  involving  Bessel  functions.  However,  their  solution  involved  a 
very  complicated  and  a tedious  computational  procedure.  They  demonstrated 
that  the  resonance  condition  similar  to  that  of  the  wind-tunnel  wall  inter- 
ference problem,  may  occur  in  staggered  cascades  with  arbitrary  interblade 
phase-lag,  the  resonance  parameter  values  depending  upon  Mach  number, 
oscillatory  frequency,  interblade  phase-lag,  gap,  stagger,  and  acoustic 
velocity. 

Fleeter^i  studied  the  effects  of  compressibility  on  both  fluctuating 
lift  and  the  fluctuating  moment  coefficients  for  cascaded  airfoils  due  to 
an  upstream  nonuniformity.  Through  the  application  of  Fourier  transform 
theory.  Fleeter  reduced,  the  time-dependent,  compressible,  two-dimensional 
partial  differential  equation  which  describes  the  perturbation  velocity  ^ 
potential,  into  an  integral  solution  equation.  Although  the  details  of  the 
solution  were  not  given  in  the  paper.  Fleeter  seems  to  have  adopted  the 
laborious  procedure  used  by  Lane  and  Friedman. 20 

Kaji  and  Okazaki^z  investigated  the  problem  of  propagation  of  sound 
waves  through  a blade  row  using  the  acceleration  potential  method.  They 
made  use  of  the  singular  point  method  in  a rigorous  and  general  form,  in 
which  the  kernel  function  of  the  aerodynamic  integral  is  expressed  by  a 


rapidly  convergent  series  of  exponential  functions. 

Rao  and  Jones^^  adopted  the  present  lifting  surface  technique  for 
determining  the  airload  and  moment  coefficients  on  a typical  airfoil  of 
a staggered  cascade  of  airfoils  in  subsonic  flow.  Circumferential 

distortion  due  to  inflow  conditions  is  expressed  as  an  interblade  phase- 
lag,  and  both  oscillatory  inflow  and  the  case  of  oscillating  blades  are 
considered.  Results  are  obtained  for  several  values  of  frequency,  Mach 
number,  and  interblade  spacing,  stagger,  and  phase-lags.  The  results 
of  the  oscillatory  flow  case  compare  well  with  those  of  Fleeter. 

Runyan  and  Woolston^*^  solved  for  the  aerodynamic  loads  on  oscillating 
wings  by  modifying  the  concepts  of  Falkner^  for  steady  lifting  surface 
theory.  The  loading  on  the  wing  was  determined  by  solving  a set  of 
simultaneous  equations  that  were  generated  by  representing  the  loading 
with  a series  representation  compatible  with  boundary  conditions.  In 
another  approach,  Lawrence  and  Cerberus  simplified  the  governing  equations 
by  using  the  assumption  of  low  aspect  ratio.  They  approximated  the  velocity 
potential  equation  by  using  a weighting  factor  in  the  integral,  and  then 
replacing  the  integrals  with  approximate  functions.  Comparison  with  the 
results  of  references  26  and  27  indicated  that  the  Lawrence  and  Gerber 
method  was  only  valid  for  wings  with  aspect  ratio  less  than  4.  The  main 
disadvantage  of  the  method  is  that  it  cannot  predict  the  pressure  dis- 
tribution. Dengler  and  Goland^s  extended  the  lifting  line  formulation 
of  Weissinger  to  the  oscillatory  case.  Their  method  is  much  simpler 
mathematically  than  lifting-surface  methods;  however,  the  chordwise 
pressure  distribution  cannot  be  predicted. 

The  pressure  distribution  measured  for  pitching  and  vertical  oscillations 
of  low  aspect  ratio  wings  by  Laidlaw  and  Halfman^s  were  compared  to  the 


results  of  Lawrence  and  Gerber^s.  This  comparison  verified  the  validity 
of  Lawrence  and  Gerber  method  for  low  aspect  ratio  wings. 

The  references  cited  above  have  adopted  the  kernel  function  repre- 
sentation and  have  illustrated  several  disadvantages.  First,  the  com- 
plexity of  the  kernel  function  involved  in  the  integral  equation 
necessitates  approximations  to  obtain  a solution.  Second,  the  kernel 
function  method  is  difficult  to  use  numerically  because  of  the  singu- 
larities and  rapid  variations  near  the  leading  edge  of  the  airfoil. 

On  the  other  hand,  the  present  lifting  surface  method  based  on  the  velocity 

potential  formulation  is  conceptually  very  simple,  and  hence  can  be 
adapted  to  any  arbitrary  planar  thin  surface(s).  Also,  in  the  case  of 
compressible  flow,  the  area  integration  can  be  replaced  by  line  integration 
with  the  aid  of  Biot-Savart's  law,  thereby  reducing  the  computational 
time  substantially.  Unlike  the  kernel  function  method,  the  present  lifting 
surface  technique  deals  with  the  circulation  which  is  more  suitable  for 
numerical  solutions. 

The  discussion  so  far  has  been  limited  to  the  application  of  lifting 
surface  methods  to  two-dimensional  and  three-dimensional  rectilinear  flows. 

What  follows  is  a brief  account  of  the  unsteady  airloads  prediction 

methods  for  helicopter  rotor  blades. 

More  papers  have  been  published  on  the  subject  of  unsteady  aerodynamics 
of  rotor  blades  during  the  past  few  years  than  in  all  the  earlier  years  of 
the  helicopter  development.  In  earlier  work,  attention  was  focused  on  vi- 
bration problems  which  often  limit  helicopter  performance.  With  the  advent 
of  high  speed  helicopters,  it  has  also  become  important  to  study  the  com- 
pressibility effects  and  the  blade  flutter  problems.  Jones  et  al^o  gave 
a detailed  account  of  significant  developments  in  the  field  of  unsteady 
aerociynamics  of  helicopter  rotor  blades.  Specifically,  the  latest  advances 
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in  such  problem  areas  as  stall  flutter  of  a retreating  rotor  blade,  flutter 
of  the  advancing  blade,  transient  effects  due  to  the  interaction  of  the  tip- 
vortex  of  one  blade  with  a following  blade  and  wake  induced  instabilities 
in  hovering  and  low-speed  flight,  were  discussed. 

Problem  areas  of  blade  loading  and  the  transient  effects  due  to  blade 
interaction  with  the  wakes  of  other  blades  and  the  flutter  characteristics 
of  the  advancing  blade  are  considered  here.  It  has  been  found  that  under 
certain  operating  conditions,  rotor  blades  can  flutter  in  both  hovering 
and  forward  flight.  This  phenomenon  has  been  investigated  by  several 
researchers  and  their  studies  are  of  particular  interest.  It  has  been 
known  for  a long  time  that  the  proximity  of  the  helical  wake  is  a con- 
tributing factor  to  blade  flutter.  When  the  rotor  has  high  inflow  through 
the  disk,  the  wake  of  a blade  is  removed  rapidly  away  and  its  effect  on 
the  aerodynamic  forces  on  the  blade  is  approximately  the  same  as  that  of  a 
rectilinear  wake.  For  the  case  of  low  inflow,  the  wake  is  closely  coiled 
under  the  rotor  disk  and  it  can  have  a strong  influence  on  the  aerodynamic 
forces  acting  on  a blade.  For  low  inflow  conditions,  Loewy^i  used  a 
simplified  mathematical  model  of  the  flow  and  developed  a general  theory. 

He  was  able  to  derive  formulae  for  the  aerodynamic  forces  on  a typical 
blade  section  as  a function  of  reduced  frequency,  the  frequency  ratio, 
and  wake  spacing.  The  special  function  that  he  obtained  was  similar  to 
that  of  Theodorsen's  function  for  a single  airfoil  in  straight  flow. 

He  has  shown  that  the  airloads  are  highly  oscillatory  for  the  smaller 
values  of  wake  spacing.  Since  the  stability  of  the  rotor  blade  depends 
on  its  aerodynamic  coefficients,  Loewy  also  investigated  the  variation  in 
the  pitching  moment  damping  coefficient  of  a particular  blade  section  as 
the  frequency  ratio,  p/n,  varied  for  specified  positions  of  axis  of 
oscillation  and  a range  of  values  of  wake  spacing.  He  found  that  the  damping 
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coefficient  became  negative  whenever  p/0  was  slightly  greater  than  an 
integer  for  axis  oscillation  either  forward  or  aft  of  quarter-chord. 

Similarly,  he  found  that  the  damping  coefficient  for  a flapping  oscil- 
lation dropped  sharply  at  integral  values  of  p/fi  but  did  not  actually 
become  negative.  From  these  results,  he  concluded  that  at  low  inflow, 
wake-excited  pitching  oscillations  of  the  blade  might  occur  when  the 
frequency  ratio  is  approximately  an  integer  and  the  pitch-axis  is  not 
at  quarter-chord. 

Whereas  Loewy  developed  his  theory  for  a multi bladed  rotor,  J.  P. 

Jones 32  had  independently  treated  the  simpler  case  of  a single  rotating 
blade  in  hovering  flight.  He  applied  a method  developed  by  W.  P.  Joneses 
to  derive  the  appropriate  aerodynamic  coefficients  for  an  oscillating  single 
rotor  blade  for  use  in  his  flutter  analysis.  He  approximated  the  actual 
flow  conditions  by  neglecting  curvature  effects  and  assuming  a simple  two- 
dimensional  mathematical  model  consisting  of  a reference  blade  and  its 
wake  and  an  infinite  number  of  wakes  lying  beneath  the  reference  blade 
extending  from  -“  to  <».  He  considered  flapping  and  pitching  motions  and 
compared  his  results  with  those  obtained  experimentally  by  Daughaday  and 
Kline 3**.  On  the  basis  of  this  work,  it  was  concluded  that  the  wake  is 
primarily  responsible  for  some  of  the  vibratory  phenomena  found  on  helicopters 
in  practice. 

Timman  and  Van  de  Vooren35,  on  the  other  hand,  assumed  that  there  was 
no  inflow  through  the  rotor  disk  and  developed  a theory  for  calculating  the 
aerodynamic  forces  on  a blade  rotating  through  its  own  wake.  Their  results 
agree  with  those  obtained  by  Loewy3i  and  Jones32  in  the  limit  when  zero 
spacing  between  the  wakes  is  assumed.  Flutter  calculations  were  done  for 
a two-bladed  rotor' model  and  their  predictions  were  compared  with  experimental 
observations.  They  found  that  wake  induced  flutter  could  occur  under 
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certain  conditions.  When  the  inertia  axis  of  the  blade  is  aft  of  its 
pitching  axis,  wake  flutter  is  possible  if  the  natural  flapping  freijuency 
is  larger  than  the  natural  pitching  frequency.  Conversely,  when  the  inert! 
axis  is  forward,  the  pitching  frequency  has  to  be  higher  than  the  flapping 
frequency  for  flutter  to  occur.  All  this  theoretical  work  confirms  the 
conclusion  that  the  proximity  of  the  wake  is  a contributing  factor  to 
rotor  blade  flutter. 

All  the  theoretical  work  described  above  is  based  on  the  assumption 
that  the  flow  is  incompressible.  However,  with  the  advent  of  helicopters 
capable  of  flying  with  blade  tip  speeds  ranging  up  to  and  in  excess  of 
the  speed  of  sound,  compressibility  effects  need  to  be  taken  into  account 
when  determining  aerodynamic  coefficients  for  use  in  flutter  analysis. 

Jones  and  Rao^^  were  able  to  do  this  and  have  computed  aerodynamic 
coefficients  for  a range  of  Mach  numbers,  reduced  frequencies,  and  wake 
spacing.  Their  analysis  is  based  on  the  application  of  the  theory  of 
Jones37,  for  a single  airfoil  oscillating  in  compressible  flow,  to  the 
Loewy  two-dimensional  mathematical  model  of  the  helical  wake.  The  values 
of  the  aerodynamic  coefficients  agree  with  those  obtained  by  Loewy 
and  J.  P.  Jones 32  for  zero  Mach  number  but  differ  appreciably  as  the  Mach 
number  is  varied. 

Hammond38  also  developed  a theory  for  determining  compressibility 
effects  by  using  a different  model  of  flow  from  that  used  in  Ref.  36.  In 
his  model,  the  wake  of  the  qth  blade  of  a Q bladed  rotor  after  n revo- 
lutions extends  from  -2Tr(n+q/Q)  to  «>;  in  Jones  and  Rao's  model  it  extends 
from  -<»  to  ».  His  aerodynamic  coefficients  for  several  Mach  number  and 
inflow  ratios  are  in  general  agreement  with  the  results  of  Jones  and  Rao. 

While  the  aerocjynamic  derivatives  predicted  by  two-dimensional  strip 
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theory  are  widely  used  in  predicting  the  flutter  speeds  of  helicopter 
rotor  blades,  the  method  does  not  allow  for  curvature  and  finite  aspect 
ratio  effects.  For  incompressible  flow,  Ashley  et  al  in  Ref.  39  developed 
a three-dimensional  model  in  which  they  modified  Reissner's  theory**®  for 
oscillating  wings  in  rectilinear  flow  by  including  the  free  stream  velocity 
variations  along  the  span.  Their  results  indicate  a negligible  difference 
between  two  and  three-dimensional  solutions  up  to  95%  of  the  span.  Jones 
and  Rao**!  also  studied  tip-vortex  effects  in  compressible  flow  and  they 
al  so  concluded  that  such  effects  are  negligible  except  in  regions  close  to 
the  tip. 

The  earlier  attempts  to  include  the  tip  vortices  in  a truly  three- 
dimensional  fashion  involved  a helical  wake  representation  below  the 
rotor  blade  in  hover  flight  and  a skewed  helical  pattern  for  forward 
flight.  Miller**^  summarized  some  of  the  work  that  has  been  done  at  M.I.T. 
and  elsewhere  up  to  that  time  on  the  prediction  of  unsteady  airloads  on 
helicopter  rotor  blades.  He  showed  that  by  replacing  the  curved  tip  vortex 
from  the  preceding  by  a straight  line  vortex  which  is  tangential  to  the 
tip  vortex  at  the  point  where  it  crosses  under  a following  blade,  he  could 
obtain  reasonable  estimations  of  its  effect.  It  was  pointed  out  that  almost 
all  the  harmonic  content  of  the  blade  air  load  of  interest  in  rotor  vibration 
is  that  contained  in  harmonics  above  the  second.  By  representing  the 
oscillatory  wake  of  the  blade  itself  in  a relatively  simple  way  and  treating 
the  tip  vortex  of  the  preceding  blade,  he  was  able  to  obtain  reasonable 
results. 

43  44 

In  some  of  his  earlier  work.  Miller  ’ developed  a helical  wake 
model  in  which  the  rotor  wake  was  divided  into  a "near"  wake  and  a "far" 
wake.  The  near  wake  included  the  portion  attached  to  the  blade  that 


extends  approximately  one-quarter  of  a revolution  from  the  blade  trailing 
edge.  The  effects  of  the  near  wake  include  an  induced  chordwise  variation 
in  downwash  and  were  formulated  using  an  adaptation  of  Loewy's  strip 
theory.  The  chordwise  variation  in  the  velocity  over  the  airfoil  induced 
by  the  far  wake  was  neglected.  Miller  extended  his  model  to  study  the 
forward  flight  case  and  found  that  the  nonuniform  downwash  induced  at  the 
rotor  disk  by  the  wake  vortex  system  could  account  for  the  higher  harmonic 
airloads  encountered  on  rotor  blades  in  forward  flight.  He  also  showed 
that  under  certain  conditions  of  low  inflow  and  low  speed  transition  flight 
the  returning  wake  could  be  sucked  up  into  the  leading  edge  of  the  rotor 
which  would  account  for  some  of  the  vibration  and  noise. 

Several  investigators  used  better  representations  of  the  helical  or 
skewed  helical  wake  for  steady  flight.  Jenney  et  al'*^  reviewed  several 
simpler  conventional  methods  for  computing  hover  performance  and  pointed 
out  that  one  of  the  primary  reasons  for  discrepancies  in  results  are  due  to 
the  use  of  inadequate  mathematical  wake  models  in  the  conventional  theories. 
They  attributed  the  major  discrepancies  to  wake  contraction  in  close 
proximity  to  the  rotor.  They  developed  a theoretical  lifting  line  Hover 
Performance  method  representing  the  wake  by  a finite  number  of  vortex 
filaments.  They  prescribed  the  wake  geometry  in  the  computer  program  in 
which  provisions  were  made  to  include  a wide  variety  of  wake  geometries 
ranging  from  the  classical  uncontracted,  momentum  axial  velocity  geometry 
(classical  wake)  to  the  more  realistic  contracted,  variable-axial-velocity 
geometry. 

Landgrebe"*^ subsequently,  conducted  an  analytical  and  experimental 
Investigation  to  acquire  systematic  model  rotor  performance  and  wake 
geometry  data  and  to  evaluate  the  accuracy  of  various  analytical  methods  in 
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predicting  the  effects  on  performance  of  changes  in  helicopter  rotor  design 
and  operating  parameters.  A rapid  computer  analysis  was  developed  for 
computing  the  distorted  wake  geometry  of  a rotor  wake.  The  numerical 
procedure  involved  two  basic  steps  - the  use  of  the  Biot-Savart  law  to 
compute  the  wake  distortion  velocities  produced  by  s given  wake  geometry 
and  the  integration  of  these  velocities  over  a small  time  increment  to 
establish  new  wake  geometry.  These  steps  were  successively  repeated  until 
a converged  periodic  wake  geometry  resulted.  It  was  found  that  analysis 
based  on  a contracted  wake  geometry  generally  provided  significantly  im- 
proved predictions  of  performance  for  those  rotor  operating  conditions  whore 
the  more  classical  uncontracted  wake  analyses  exhibited  major  shortcomings. 
Through  extensive  correlation  testing,  he  validated  the  general  trends 
predicted  by  his  numerical  work.  Analysis  of  the  wake  data  for  that  portion 
of  the  wake  which  was  stable  (i.e.  near  the  rotor)  indicated  that  the 
data  could  be  expressed  in  relatively  simple  generalized  equations  which 
facilitate  the  rapid  estimation  of  contracted  wake  geometries  for  a wide 
range  of  rotor  design  and  operating  conditions. 

In  an  independent  study.  Piziali**®,  developed  an  alternative  numerical 
method  in  which  the  wake  of  the  rotor  blade  is  represented  by  discrete 
straight  lines  and  trailing  vortex  elements.  He  satisfied  the  chordwise 
boundary  conditions,  but  the  rotor  blade  was  limited  to  one  degree-of- 
freedom  in  flapping.  His  method  is  also  limited  by  the  use  of  empirical 
factors  which  adjusts  the  discrete  time  positioning  of  the  shed  vortices. 

Sadler**®,  using  a model  similar  to  that  of  Piziali,  developed  a 
method  for  predicting  the  helicopter. wake  geometry  at  a "start  up"  config- 
uration. He  represented  the  wake  by  a fine  mesh  of  transverse  and  trailing 
vortices  starting  with  the  first  movement  of  the  rotor  blade  generating  a 
bound  vortex,  and,  to  preserve  zero  total  vorticity,  a corresponding  shed 
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vortex  in  the  wake.  Integrating  the  mutual  inteference  of  the  trailing 
and  shed  vortices  upon  each  other  over  small  intervals  of  time,  Sadler 
was  able  to  predict  a wake  geometry.  Although  his  model  showed  fair 
agreement  with  the  available  experimental  data  for  advance  ratios  above 
one-tenth,  Sadler's  method  is  limited  due  to  the  large  computational  time 
required  to  represent  the  wake  by  a finite  mesh. 

Another  important  problem  is  the  prediction  of  the  unsteady  airloads 
in  normal  forward  flight  and  allowing  for  blade-tip  vortex  interaction. 
Landgrebe  and  other  investigators  demonstrated  that  as  a given  blade 
advances  the  tip  vortex  from  the  preceding  blade  may  pass  closely  under- 
neath it  and  the  point  of  intersection  moves  inward  along  the  blade  as  'i 
increases.  Ward  and  Snyder^®  predicted  the  extent  of  the  inward  movement 
as  a function  of  the  advance  ratio  and  the  number  of  blades.  In  their 
study,  they  did  not  account  for  the  actual  induced  response  due  to  the 
movement  of  the  tip  vortex  but  instead  considered  the  analogous  problem  of 
the  response  of  the  blade  to  a repetitive  moving  load.  Their  results 
showed  for  flight  at  high  advance  ratios  (ii>0.25)  the  main  effect  of  the 
moving  load  was  to  amplify  response  in  the  lower  resonance  modes. 

The  topic  of  the  evaluation  of  appropriate  aerodynamic  coefficients 
for  use  in  the  analysis  of  classical  bending- torsi  on  flutter  of  rotor  blades 
is  of  great  interest  to  a helicopter  designer  in  view  of  the  present  day 
helicopters  with  high  tip  speeds.  Shipman  and  Wood^i  considered  this 
problem  but  they  did  not  take  compressibility  and  finite  aspect  ratio  ef- 
fects into  account.  The  two-dimensional  mathematical  model  used  was 
similar  to  that  employed  by  other  authors  except  that  they  assumed  that 
flutter  would  first  occur  when  the  relative  velocity  over  the  rotor  blade 
reaches  its  critical  value  when  'r  = tt/2.  For  greater  or  lower  values  of 
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the  realttve  speed  would  be  reduced  below  the  critical  speed  for  flutter 
and  any  incipient  growing  flutter  oscillation  would  be  damped.  This 
reasoning  led  them  to  represent  the  blade  motion  by  a series  of  oscillatory 
pulses,  where  each  disturbance  occurs  over  the  range,  Tr/2  - A4'^<  f < ir/2  - ai'2- 
Corresponding  to  each  burst  of  oscillation,  packets  of  vorticity  are  assumed 
to  be  shed  in  the  wake.  With  increasing  forward  speed,  the  spacing  between 
the  packets  of  vorticity  also  increases  and  it  was  found  that  the  flutter 
speed  became  constant  when  p,  the  advance  ratio,  was  above  0.2. 

The  approach  used  in  a study  by  Rao  and  Jones^^  differs  from  that 
adopted  by  Shipman  and  Wood  in  that  continuous  high  frequency  small 
oscillations  are  assumed  to  be  superimposed  on  the  normal  periodic  motion 
of  the  blade.  Rao  and  Jones  developed  a simple  but  general  numerical 
lifting  surface  method  for  predicting  unsteady  airloads  on  a single-bladed 
rotor  blade  on  a full  three-dimensional  basis.  The  numerical  method  was 
based  on  the  velocity  potential  formulation  and  was  not  restricted  in  any 
way  as  to  frequency,  mode  of  oscillation  or  aspect  ratio  when  M < 1. 

However,  in  their  model  they  used  a classical  wake  model  and  considered 
uniform  downwash  only.  The  general  theory  was  developed  for  a rotor  blade 
at  the  v = Tr/2  position  where  flutter  is  most  likely  to  occur.  Calculations 
of  aerodynamic  coefficients  for  use  in  flutter  analysis  were  made  for 
forward  and  hovering  flight  with  low  inflow  for  Mach  numbers  0 and  0.8 
and  frequency  ratios  p/n  = 1 and  4.  The  results  were  compared  with  values 
given  by  two-dimensional  strip  theory  for  a rigid  rotor  hinged  at  its  root. 

The  comparisons  indicated  the  inadequacies  of  strip  theory  for  airload 
prediction.  One  important  conclusion  drawn  from  this  study  was  that  the 
curved  wake  has  a substantial  effect  on  the  chordwise  load  distribution. 

The  pitching  moment  aerodynamic  coefficients  differed  appreciably  from  the 
results  given  by  strip  theory. 
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II.  AERODYNAMIC  THEORY 

Governing  Equation  of  Flow 

Consider  an  airfoil  in  motion  in  a continuous  fluid  medium.  A 

I Cartesian  coordinate  system  is  chosen  such  that  the  origin  0 lies  at 

[ 

I 

I the  midchord  point  of  the  airfoil,  and  the  positive  x-axis  in  the 

direction  of  the  free  stream  relative  to  the  airfoil.  The  system  of 
coordinates  is  shown  in  the  sketch  below. 


Velocity  of  a fluid  particle  at  any  point  P(Xp,  yp,  Zp)  in  the  flow 


field  is 


“V  = “l.,  t + <V,  J + k . ( 

In  a potential  flow  around  the  airfoil,  the  laws  of  conservation  of 
mass  and  momentum  are,  respectively. 


(2.1) 


- + =0  , 


(2.2) 


P K 4 

— + 1 vb  = 0 

pt  ^ 

Expressing  the  velocity  in  terms  of  a potential  function  , 
the  continuity  equation  (2.2)  reduces  to 

^ /»  p-t  * 

For  a barotropic  fluid,  p = p(p)  and,  hence 

rf  _ __  _j Pb 

* /»a^  pt 

Along  a streamline,  the  momentum  equation  yields 

El  ^ ^ _ JL  Ek.  - I Dt* 

vt'-z'  /’9t  "/’pt  ■ 

Also,  from  the  momentum  equation, 

vF  l-ji)  + 1 = 0 


(2.3) 


(2.4) 


(2.5) 


(2.6) 


(2.7) 


Equations  (2.5),  (2.6)  and  (2.7)  are  combined  to  obtain  the  governing 


L 
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potential  flow  equation, 


(2.8) 


This  general  form  of  the  equation  is  applicable  to  unsteady  compressible 
flows  with  no  shock  waves.  Based  on  the  type  of  flow,  equation  (2.8) 
is  modified  as  follows. 

a.  Incompressible  Flow  (Steady  and  Unsteady  State)- 

Local  speed  of  sound  in  an  incompressible  medium  is  infinite  and 
he.ice  the  governing  differential  equation  of  potential  flow  becomes 
Laplace's  equation, 


^ = 0 . (2.9) 

b.  Compressible  Flow  (Unsteady  State)- 

Velocity  of  a fluid  particle  in  the  vicinity  of  the  airfoil  can  be 
written  as 

+ vrj  , (2.10) 


where  the  perturbation  velocity  components  u,  v,  and  w are  considered 
small  in  comparison  to  the  free  stream  velocity  component.  Coordinate 
transformations 

x=xi 
> = Yl/(i 
Z = Z£/P 
t = T£/U 

P - , 


(2.11) 


where. 
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and  the  substitution 

If  = VI  f CX.Y.Z)  e ’2) 

modify  equation  (2.8)  to 


V*  $ + 

= 0 . 

(2.13) 

' - ^ 

i- 

(2.14) 

K = Mo/p’’  , 

(2.15) 

5k  = Mk  . 

(2.16) 

c.  Compressible  Flow  (Steady  State). 

Frequency  of  external  excitations  p'=0  in  the  steady  state  resulting 
In  100,  and  the  governing  flow  equation  becomes  Laplace's  equation  in  the 
transformed  coordinates. 


(2.17) 


Velocity  Potential  in  Terms  of  Boundary  Values 

To  obtain  a relation  for  the  velocity  potential  $ in  terms  of  its 
value  and  the  normal  derivative  on  known  boundaries  so  as  to  formulate  a 
boundary  value  problem  in  conjunction  with  Helmholtz's  equation  (2.13),  we 
now  refer  to  the  sketch. 


Fluid  region  D surrounding  the  airfoil  is  bounded  externally  by 
the  surface  z,  and  encloses  an  impermeable  surface  S comprising  the 
airfoil  and  the  attached  wake  surfaces. 

With  the  use  of  Gauss'  divergence  theorem  for  the  vector  B , 
a reciprocal  theorem  can  be  proved  as 


^ f |i  js  = 0 |i  as  . 

where  $ and  both  satisfy  Helmholtz's  equation  (2.13),  and  n represents 


(2.18) 


the  unit  normal  vector  at  any  point  on  S and  z,  directed  toward  the  fluid. 


where  R is  the  distance 


Assuming  where  R is  the  distance 

between  a flow  point  T’Cx^.y|..zp  and  a point  <5  cx  :y,z) 
on  (S  + e).  In  order  to  avoid  the  singularity  in  i^when  the  points  P 
and  Q coincide,  P is  enclosed  in  a small  spherical  surface  a which,  in  the 
limit,  is  taken  to  the  point  P.  Further,  as  the  exterior  boundary  is 
extended  toward  infinity,  contributions  over  e in  equation  (2.18)  approach 
zero  for  point  P in  the  vicinity  of  the  airfoil.  Equation  (2.18)  is  now 
rewritten  as 

-iKR 

S -XKR 


(2.19.' 


where  X,  Y,  Z are  the  independent  variables  in  the  integrand  on  the  right 
side.  This  equation  provides  a means  to  express  the  velocity  potential 
at  any  point  P in  the  flow  in  terms  of  its  value  $ ^nd 

the  normal  derivative  M on  the  surface  of  the  airfoil  and  its 

wake.  With  the  assumption  of  thin  airfoil,  the  closed  surface  of  inte- 
gration can  be  divided  into  two  equal  and  parallel  surfaces.  Also,  the 
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(2.21) 


physical  boundary  condition  on  these  surfaces  requires  the  relative  flow 
to  follow  the  airfoil  and  its  wake  geometries.  That  is, 

SI  ( ^ J + Vi  = O ^ (2.20) 

whereby 

di 

— = -n„ 

* 

It  is  further  noted  that  on  the  upper  surface  of  integration 
a/ a -3/32  and  on  the  lower  side 

Hence,  from  equation  (2.19),  the  downwash  induced  at  a point  P on  the 
airfoil  in  a direction  along  the  positive  z-axis  is  given  by 

* x —IK  R 

47.  = 4XW^=  ]J  di  f ) ^5.  ,3. 

Qirfoii  + 
walc« 


22) 


Rewriting, 


4;cV^  = If  K(x,y,z)g5( 

wak« 


-iKR 

R > 


(2.23) 


where  K = 


defines  the  strength  of  the  doublet  distribution 
over  the  airfoil  and  wake  surfaces. 

Equation  (2.23)  is  the  most  general  form  of  the  velocity  potential 
formulation  applicable  to  unsteady  flow  of  a compressible  fluid  past  a 
thin  airfoil  at  small  angles  of  attack.  It  can  be  suitably  modified  for 
incompressible  flow  (steady  or  unsteady),  two-dimensional  compressible  flow 
(steady  and  unsteady)  and  three-dimensional  steady  compressible  flow. 

For  all  steady  cases,  and  unsteady  incompressible  flow,  tC=  0,  and 
equation  (2.23)  becomes 


47^ 


air  fall  4-  Y 


)dS  . 


(2.24) 
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For  two-dimensional  subsonic  unsteady  flow,  the  velocity  potential 
formulation  reduces  to 

aicWj,  = 11  <(x,2)  ^ , 


-I  0 

where  P is  a point  on  the  Y=0  axis,  and 


R = [(x-xp%  Y^+ 

Equation  (2.25)  is  rearranged  as 


= I I 

^ -I  *•  9Zj,  o 


-4.1CR 
S 


We  now  let 


whereby 


or 


Therefore, 


Y*= 


zydV  = 2o<*T  d'C 


*0  — ^ iC  R _ 

( I-  ay  = ( 

o rs 


dT 


where 


v/t^-1 

|.|-  [CX-xp*  + Zj!]*  .nd 


} 


(2.25) 


(2.26) 


dy]dx 

• 

(2.27) 

dnci  j. 

(2.28) 

(2.29) 

• 

(2.30) 

(2.31) 


Ho  (KMI)  = Jft(ICl-fl)  - 4.  (fC|o<|)  ^ is 

the  modified  Hankel  function  of  the  second  kind  Use  of  result  (2.31)  in 
equation  2.27  yields 

ax V,,  = - >^  J k (x.z)  Ip,  [ wf'^c K ki;]  ‘♦X  . 

2. 


(2.32) 


(2.33) 


r 
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The  function  , defined  by  equation  (2.32)  satisfies 

the  Bessel's  equation 

A t + »i%y  =0 , (2.34) 

ar  * 

where 

. (2.35) 

CO 

Hence  it  can  be  shown  that  H,  CKNI)  also  satisfies  the 
two-dimensional  wave  equation. 


^ 3^  2^  .-a,  ^ I .(-O 

^ 3zf 


+ «*  ) i«i)  =0 


(2.36) 


To  show  this,  we  consider  the  following. 

XO  . /a.) 


where 


« 

ee 


Also , 


Mo  = H.  cr^  ,(r=M^NO 
r=  . 


ax 


o = 


^ - 1.  r it  sj*  1 

~ ‘'I* 


JliL 

J|* 

<4ie’ 


and 


dX 


I* 


^ U.  /3r  ^ ^ dUfi.iii 

WX,,''  Jf*  JX^ 


it  = - 


3X 


► 


(2.37) 


(2.38) 


(2.39) 


(2.40) 


Therefore , 


^ r jc^Cx-xp)  K* 

JXf  1 — / ^ jr 

Substituting  equations  (2.40)  and  (2.41)  in  equation  (2.39), 

(j.)  r K*  _ K^cx-xt)’'! 

L r p j 


+ '•Hf 


-It* 


Similarly, 


!-[»>.]  - , 


az 


and 


I’ 


itl 

3*7  J 


Now, 


= 2!t 

Y ''  »z/  -If 


3Z,.  " I* 


^vid 


_ K.’^Zk  / K^ZbX 

r T ( TV 
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(2.41) 


(2.42) 


(2.43) 


(2.44) 


(2.45) 


Therefore 


t.  f u‘*Vu,l  = r K*Zl  1 ^ d (c*-  ^4  yz 

?z^'-  dr*  J dF'-F 

Substituting  equations  (2.42)  ancj  (2.46)  in  equation  (2.36),  we  have. 


^ 3X.*  3Zt 

= [ e]  ^ 


d r 2 K*  K^«<* 


dr’ 


r ^ 

df*  ^ f* 


+ k‘h<*> 

o 


= K^[ 


d\‘*^ 

dr* 


^ j_  dH. 


(O 


H^] 


dr 


s 0 . in  view  of  the  equation  (2.34), 


Hence, 


whereby  the  formulation  (2.33)  becomes 


«o 


^ ^"]nr(«ix-x^i;dx. 


(2.46) 


(2.47) 


(2.48) 


27 


as  Zp,  Z 0. 

A slight  modification  to  the  above  form  can  be  achieved  by  writing 


from  Bessel's  equation  (2.34).  Also  as  Zp,  Z ->•  0, 

h*  = fc*Cx-xp^  , 

and 

or 


Therefore , 


dP.  = _ dx, 

|x-x^l  ^ 


dH 


. d2i(-  _ i*li  r_  ix-X|.i  1 

dXf  df*  JX^  L )c(x-xp'’ 


^ = i fiil®  dXb-| 
jtY'-  dr  ^ dX|,  • dp  J 

= J-  d*H‘*" 

K*  dXi  ■ 

Hence,  equation  (2,49)  yields 


dXt.*  “ i^KCx-xp’ 


dX 


or 


I- 


(X'Xj,)  dXj, 


(2.49) 


(2.50) 


(2.51) 


(2.52) 
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where 


R=  [(x-xp’'+ 


Rearranging  equation  (2.55), 


0 


Consider 


^ r-)  - -II* 


(2.54) 


Equation  (2.48)  thus  changes  to 

^ ’(  K I X- I ; JX  (2  53) 

= -!^r^^H^"^CKlX-X,,|)dX  (3.54) 

“I  (X ^1?/ 

Equations  (2.53)  and  (2.54)  are  two  convenient  forms  of  the  two- 
dimensional,  unsteady,  compressible  flow  formulation. 

For  two-dimensional  incompressible  flow,  both  steady  and  unsteady 
states,  we  now  refer  to  equation  (2.24)  and  rewrite 

= J j-  K(x.7;^ /i)dy4X  X (3.55) 

-1  A cr/i.  ' 


(2.56) 


(2.57) 


f ^ 


" t [(x-xpS4"][(x-xp%Y%Z^"]"'| 


(x-x^)V 


(2.58) 
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Therefore,  equation  (2.56)  becomes# 


.»u  9 


(2.59) 


which  reduces  to 


- P<x,  ^ 


(2.60) 


as  Zp,  Z ->  0,  and  noting  at 


r 

U-  a 

r ^-4  1 

•-  (x-x^)V 

■'  ax  i 

Kx-xpV  2^  J 

A d 


Alternatively,  equation  (2.60)  is 

X-Xf, 

Three-dimensional  incompressible  flow,  steady  and  unsteady,  formulation 
is  simplified  further  as  follows.  From  equation  (2.24)  , 


(2.63) 


where 


airfoil 

wak€ 


e = [(x-x^)\  (Y-yp% 

> ^1* 

= s_  r 2l=-i^n  +i-r 

5X  L r3  J ay  L ^3  J 


(2.64) 


. (2.65) 
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Therefore,  as  Z 0,  and  for  a thin  lifting  surface  at  small  angles  of 
attach  K(Zp  ^ 0) 


WAke 


Numerical  Scheme  to  Solve  Integral  Equations 

The  integrands  of  the  surface  integrals  appearing  in  the  velocity 
potential  formulations  for  various  types  of  flows  are  functions  of  the 
unknown  discontinuity  in  the  modified  velocity  potential,  K(X,Y),  and 
known  functions  of  the  airfoil  and  wake  geometry  in  the  transformed 
coordinates.  These  also  depend  upon  the  angular  frequency  of  external 
excitations  that  are  responsible  for  the  unsteady  state  of  the  relative 
flow  past  the  airfoil.  In  order  to  facilitate  evaluation  of  K-distribution, 
the  lifting  surface  is  divided  into  a riumb’er  of  conveniently  shaped  boxes 
and  the  value  of  K over  each  box  is  assumed  to  be  constant.  Making  use 
of  the  wake  boundary  condition  that  there  can  exist  no  pressure  dis- 
continuity in  the  wake,  the  value  of  K at  any  spanwise  location  in  the 
wake  can  be  expressed  as  a function  of  the  value  of  K at  the  trailing 
edge  box  at  the  same  spanwise  location.  Hence,  the  wake  is  divided  into 
a number  of  trailing  strips  and  K(X,Y)  is  considered  constant  across 
each  strip  but  to  vary  with  X downstream.  When  the  downwash  contributions 
due  to  all  the  boxes  and  strips  are  obtained,  the  problem,  in  general, 
is  reduced  to  solving  simultaneously  a set  of  as  many  linear  equations  as 
the  number  of  boxes  on  the  lifting  surface.  Steady  and  unsteady  load 
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calculations  are  then  carried  out  in  terms  of  the  known  doublet  distribution 
and  the  lifting  surface  geometry. 


III.  FIXED  WING  APPLICATIONS 
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Two-Dimensional  Incompressible  Flow 
a.  Steady  State. 

The  relation  for  the  downwash  at  a point  on  an  airfoil  is  expressed 
as  equation  (2.60), 


00 


= f K I-  [ 


X-Xi 


-«  (X-  t 

Z-»o  ' 


] 


dX 


1) 


N 


_ Kte  + i K r * ! 1 

“ X.'i  j7i  ^'-X.-Xj^-J  x.-x.-pJ- 


The  airfoil  Is  divided  Into  N equal  strips  of  width  2D  and  K is  assumed  to 
be  a constant  over  each  strip.  is  related  to  Kj^  and  is  normally  taken 
as  equal  to  Kj^.  Equation  (3.1)  can  be  written  as 


= [a]{K}  , 


(3.2) 


where 


and 


"•J 


1 


1 


X.-X:  + P 

■t  J 


x.-x. -P 

^ J 


1 


X.  -1  +zp 

A0 


I 
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For  a prescribed  boundary  condition,  is  known.  This  is  a non- 
dimensional  downwash  velocity  at  a point  X^.  . For  a constant  angle  of 
attack  case,  = a(i  = 1,  2,  ...,  N)  and  the  lift  per  unit  span  is 

l'  = />Uk.  = = />vUk 


te  - ' - N 


or 


= lVc/’u*x;  = . (3.3) 

Pitching  moment  per  unit  span  about  the  mid-chord  axis  is  (nose-up  positive) 

I 


m'= 

-I 


or 


c,.,  = M 7c /» u V;  = K..-  2 


N 


_p  z: 

H rt=|  ” 


(3.4) 


The  analytical  solution  for  K distribution  is  given  by 

K = 2.  0 -t*  siw  6)  0< 

where  e=  e.as~'C-X) 


(3.5) 


The  computer  program  based  upon  the  constant  K distribution  over 
each  strip  is  very  simple,  especially  when  the  collocation  points  of  each 
strip  are  taken  at  the  center  of  each  strip.  For  this  case,  the  results 
are  presented  below  for  several  values  of  N and  are  compared  with  the 
exact  results. 


Exact 

Lifting 

Surface  Solutions 

Sol  uti  on 

N=5 

N=20 

N=30 

6.2832 

6.2831 

6.2829 

6.2818 

^ma 

-3.1416 

-2.6431 

-2.9885 

-3.0384 
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The  computational  time  on  AMDAHL  470  for  all  three  cases  (N=5,  20, 
and  30)  was  less  than  0.54  seconds  using  WATFIV  compiler  which  is  only 
1/5  th  as  fast  as  FORTRAN  G or  H compiler.  As  can  be  seen,  the  results 
for  the  coefficients  are  in  excellent  agreement  indicating  that  the  com- 
puted K|^  is  close  to  the  exact  value  even  for  N=5  case.  The  moment 
coefficient  continues  to  reach  the  exact  value  as  N is  increased.  The 
K distribution  is  plotted  in  Fig.  i and  it  can  be  seen,  it  is  substan- 
tially different  from  the  exact  distribution,  especially  near  the  leading 
edge  region.  To  improve  this  situation  the  collocation  point  of  the 
leading  edge  strip  is  varied  from  ID  to  0.75D  from  the  leading  edge  and 
the  results  of  K(K/2Troi)  are  presented  below  for  a 5 strip  model. 

K/2Tra 


x = 

= -0.8 

-0.4 

0 

0.4 

0.8 

Exact 

0.396 

0.661 

0.818 

0.923 

0.986 

l.OOD 

0.492 

0.711 

0.852 

0.945 

1.000 

0.95D 

0.468 

0.697  ■ 

0.840 

0.935 

0.990 

0.90D 

0.443 

0.684 

0.829 

0.925 

0.980 

0.85D 

0.418 

0.670 

0.818 

0.914 

0.970 

0.80D 

0.394 

0.656 

0.806 

0.904 

0.960 

0.75D 

0.369 

0.643 

0.795 

0.893 

0.950 

As  can  be  seen  from  the  above  Table 

, the  K distribution  has  improved 

towards  the  leading  edge  as 

the  collocation  point 

location 

is  moved  toward 

the  leading  edge; 

however. 

the  trailing 

edge  values 

are  in 

error  for  this 

collocation  point 

location , 

resulting  in 

a decrease 

in  the 

lift  coefficient 

Another  scheme  is  tried  to  correct  this  situation.  In  addition  to  the 
leading  edge  strip  collocation  point  variation, the  trailing  edge  value. 
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K^g,  is  expressed  by  the  relation 


6 


in  place  of  K^g  = Kj^,  and  the  results  are  tabulated  below. 


K/2Tra 


X 

= -0.8 

-0.4 

0 

0.4 

0.8 

Exact 

0.396 

0.661 

0.818 

0.923 

0.986 

l.OOD 

0.499 

0.722 

0.866 

0.964 

1.023 

0.96D 

0.480 

0.711 

0.857 

0.955 

1.015 

0.92D 

0.460 

0.700 

0.843 

0.947 

1 .007 

0.88D 

0.440 

0.689 

0.839 

0.939 

0.998 

0.84D 

0.420 

0.678 

0.830 

0.931 

0.990 

0.80D 

0.400 

0.667 

0.821 

0.922 

0.982 

As  can  be  seen  from  the  above  table  for  the  collocation  point  lo- 
cation at  0.800D  and  the  modified  trailing. edge  relation  resulted  in  an 
excellent  agreement  with  the  exact  result.  In  a recent  paper,  Kocurek  and 
Tangier®^  discussed  a systematic  means  of  the  collocation  point  selection. 


b.  Unsteady  State. 

In  the  case  of  unsteady  flow,  a periodic  motion  is  assumed  and  the 
doublet  distribution  is  assumed  to  be  K^ex|?C  = K • 

The  relation  for  the  downwash  (amplitude)  at  a point  on  an  airfoil  is 

expressed  as  equation  (2.62) 

a„W,'= 


6) 


w 
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The  boundary  condition  in  the  wake  and  at  the  trailing  edge  is  expressed 
as 


9K‘ 

3X 


==  O 


1 


(3.7) 


and  K'  in  the  wake  can  be  expressed  in  terms  of  by  solving  equation  (3.7), 


1^/  -i.w(X-0 


(3.8) 


The  airfoil  is  divided  into  N equal  strips  and  K is  assumed  to  be  a 
constant  over  each  strip.  Equations  (3.6)  and  (3.8)  are  combined  to  yield 


00  -X4>(X-0 


dX 


X;-f 


I 

+ ^ It'.  r__! « 1 

jsi  J ^ — Xj  +•  J)  X^*  — Xj  — P J 

Making  use  of  the  wake  boundary  condition  at  the  trailing  edge  several 
possible  relations  can  be  obtained  for  in  terms  of  Kj^.  One  such 

frequently  used  relation  is 


(3.9) 


where 


= F Xn  , 

Fs  iJikoVj 


(3.10) 


Equations  (3.9)  and  (3.10)  are  combined  to  yield 


2a{w}  = [aKk} 


where 


I 


a.,  = 


I 


Xj  - Xj  - P 


/ j + N 


I 


I 


X - Xj  4-  P 


Xj-I 


X,-Xj-P 

I X-  Xi  * 


and 
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The  infinite  integrals  in  the  term  are  evaluated  by  separating 

into  real  and  imaginary  parts.  These  terms  are  expressed  in  terms  of 
infinite  series  and  are  analytically  integrated.  The  relation  is  given 
by 


I 


X-Xi 


ioO-X;)  r . n / ^ 

= fi  l{-r-ftvit^O-Xi) 

2.2j  ^ ' J 


(3.12) 


The  above  developed  formulation  is  equally  valid  for  oscillating 
airfoils  as  is  for  oscillatory  flows  (like  gust).  The  boundary  con- 
ditions are  usually  prescribed  in  either  case.  For  airfoils  undergoing 
flapping  and  pitching  oscillations  (about  the  mid-chord)  the  boundary 
condition  is 


and  for  oscillatory  flow 

, , ^ioyXi 

v/  = 


(3.13) 


(3.14) 


where  Z'  and  a'  are  the  flapping  and  hitching  amplitudes,  respectively. 

The  unsteady  lift  and  moment  coefficients  are  then  expressed  as 

/ 


= c- z + c„  •< 


and 


M 


- SS  C 1 c 0( 


(3.15) 
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where 


and 


*C  ) ^eC  ? j*i  J(.;  ' 

N 


H 


V (3.16) 


The  parentheses,  ( ),  in  equation  (3.16)  can  be  replaced  either  by  Z or  a 
and  the  appropriate  boundary  conditions  are  applied. 

The  analytical  expressions  for  oscillating  airfoils  are 

Sz=  , 

C„J=  i*oxc 


and 


c„,,=  x[c(i4 


and  for  the  gust  flow, 


where 


> 


c = 


j,-  4,7; 


h/  0^) 


and 


H 


C^'> 


(3.17) 


(3.18) 


0 - -6 to 

For  collocation  points  located  at  the  centers  of  each  strip  results 
are  presented  below  for  both  cases  along  with  the  comparison  with  the 
exact  results 


Oscillating  Airfoils 
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u 

N 

^mz 

Cm 

ma 

0.1 

exact 

0.077+0. 523i 

5.281-0.507i 

0.054+0. 261 i 

2.645-0.568i 

0.1 

10 

0. 074+0. 525i 

5.296-0. 500i 

0.056+0. 2 79 i 

2.828-0. 588i 

0.1 

20 

0. 074+0. 524i 

5.293-0.493i 

0. 055+0. 272i 

2.750-0.582i 

0.1 

30 

0. 075+0. 524i 

5.291-0.493i 

0. 055+0. 269i 

2.718-0.579i 

0.2 

exact 

0. 111+0. 914i 

4.690-0. lOOi 

0.11 9+0. 457i 

2.361-0.678i 

0.2 

10 

0.1 03+0. 922i 

4.718-0.102i 

0. 122+0. 490i 

2.524-0.689i 

0.2 

20 

0.1 04+0. 920i 

4.713-0.085i 

0. 121+0. 476i 

2.454-0.692i 

0.2 

30 

0.1 05+0. 91 9i 

4.710-0.083i 

0. 121+0. 470i 

2.426-0. 690i 

0.3 

exact 

0. 055+1. 253i 

4.347+0.443i 

0. 169+0. 627i 

2. 209-0. 721 i 

0.3 

10 

0. 043+1. 269i 

4.386+0.428i 

0.1 70+0. 672i 

2. 360-0. 71 8i 

0.3 

20 

0. 044+1. 265i 

4.378+0.456i 

0. 172+0. 652i 

2.294-0. 730i 

0.3 

30 

0. 045+1. 263i 

4.374+0.459i 

0.1 72+0. 644i 

2.269-0.731i 

0.4 

exact 

-0. 088+1. 571i 

4.134+1.005i 

0. 207+0. 785i 

2.1 30-0. 754i 

0.4 

10 

-0. 105+1. 595i 

4.1 83+0. 979i 

0.204+0. 842 i 

2.271-0.735i 

0.4 

20 

-0. 104+1. 590i 

4.173+1.015i 

0.21 0+0. 81 6i 

2.210-0. 757i 

0.4 

30 

-0. 102+1. 587i 

4.167+1.022i 

0. 211+0. 807i 

2.186-0.761i 

0.5 

exact 

-0. 312+1. 879i 

3.994+1.563i 

0. 237+0. 939i 

2.095-0.789i 

0.5 

10 

-0.332+1 .912i 

4.053+1.525i 

0.225+1. 006 i 

2.230-0.753i 

0.5 

20 

-0. 333+1. 905i 

4.039+1.570i 

0. 237+0. 975i 

2.171-0.785i 

0.5 

30 

-0.330+1.900i 

4.031+1.579i 

0.2 39+0. 964i 

2.148-0.793i 

Evaluation  of  the  above  listed  results  indicate  that  the  numerical 
solution  is  reasonably  accurate  even  for  10  strip  case.  The  computer 
program  is  very  efficient  and  the  computational  time  for  all  15  cases 
(5  id's  and  3 N's)  was  4.58  seconds  using  WATFIV  compiler.  For  oscillatory 
flow  the  numerical  results  along  with  the  comparison  with  the  analytical 
results  are  illustrated  below. 
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Oscillatory  Flow 


0) 

N 

^Aa/2n 

1 

exact 

0.3687+0.12591 

1 

5 

0.3868+0.15531 

1 

9 

0.3838+0.14571 

1 

15 

0.3812+0.13941 

1 

25 

0.3790+0.13491 

2 

exact 

0.0816+0.26801 

2 

5 

0.0707+0.28401 

2 

9 

0.0789+0.27121 

2 

15 

0.0823+0.26771 

2 

25 

0.0840+0.26641 

3 

exact 

-0.1452+0.17781 

3 

5 

-0.1705+0.15871 

3 

9 

-0.1623+0.15131 

3 

15 

-0.1562+0.15591 

3 

25 

-0.1521+0.16131 

4 

exact 

-0.1980-0.02071 

4 

5 

-0.1932-0.09481 

4 

9 

-0.2053-0.06551 

4 

15 

-0.2062-0.05091 

4 

25 

-0.2051-0.04101 

5 

exact 

-0.0812-0.15861 

5 

5 

-0.0415-0.27831 

5 

9 

-0.0645-0.19191 

5 

15 

-0.0746-0.17251 

5 


25 


-0.0794-0.16501 


The  values  presented  in  the  above  Table  are  for  very  high  reduced 
frequencies  which  are  of  considerable  importance  in  acoustics.  Also 
based  on  the  results,  one  can  conclude  that  the  numerical  lifting  surface 
technique  yields  a reasonable  solution  if  at  least  15  strips  are  considered. 
Several  other  schemes  were  attempted  varying  the  leading  edge  strip  col- 
location point  location.  These  schemes  improved  the  results  but  however, 
the  collocation  point  location  was  found  to  be  a strong  function  of 
reduced  frequency  and  this  also  increased  the  computational  time  due 
to  the  lack  of  symmetry. 
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Two-Dimensional  Coinpressible  Flow 

The  relation  between  the  downwash  at  any  point  X^.  on  a thin  airfoil 
in  subsonic  flow  and  the  doublet  distribution  is  derived  in  Chapter  II, 
equation  (2.48)  It  is  expressed  as 

W [ k"]  Hf  (k  1 x-x^O'JX.  ,3. 

The  wake  boundary  condition  is 

• ^ -n 


19) 


x^1  , 


(3.20) 


from  which  K^in  the  wake  can  be  expressed  in  terms  of  The  relation  is 


/ -Xx^CX— l) 


ia 


, x^i  . 


(3.21) 


Equations  (3.19)  and  (3.21)  are  combined  to  yield 

HX)  = X x=  lc|x-Xi) . 

The  airfoil  is  divided  into  N equal  strips  of  width,  2D,  and  K is 
assumed  to  be  a constant  over  each  strip.  Equation  (3.22)  is  then 


(3.22) 


-i-  4^-  (Xj ) + y (Xj  -p))  j-  - ^ |_ 

4-  I (X)  Jx‘}]  ,3 


23) 


I 

L 
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The  wake  boundary  condition  is  used  to  express  in  terms  of  K|^ 
and  the  relation  is 


(3.24) 


Equations  (3.23)  and  (3.24)  are  combined  to  yield  a set  of  N simultaneous 
equations,  i=l,  2,  ....  N.  The  Hankel  function  has  a singularity  when 
its  argument  is  zero.  The  evaluation  of  the  derivative  of  presents  no 
problem  since  the  point  of  evaluation  is,  at  the  least,  separated  from  the 
collocation  point  by  D.  The  derivative  of  the  Hankel  function  is 


OX 


1C 


(Xx-X) 

Ix-Xil 


Hf’'^(^lC  1 X-X,|) 


(3.25) 


The  second  term  in  equation  (3.23)  is  obtained  by  using  Simpson's 
1/3  rule  and  it  poses  a problem  since  the  ^(Xjlterm  becomes  infinite 
when  Xj  is  at  the  collocation  point,  X^.  This  situation  is  alleviated 
by  using  the  Struve's  function.  At  the  singularity  point,  the  portion 
of  the  surface  being  integrated  over  is 

KP 

\ ^Cx->X-Xjl)dX=  jr  I 

Xj-J> 


where 


and 


26 


45 


The  wake  integration  involving  an  infinite  integration  is  simplified 


and  is  expressed  in  the  form 


i-Xc  , 

6 


(3.27) 


Equations  (3.24)  - (3.27)  are  incorporated  in  equation  (3.23)  and  the  result  is 
that  the  problem  reduces  to  a set  of  simultaneous  equations  which  can  be 
solved  for  a given  set  of  boundary  conditions. 


Boundary  Conditions 

dZ  " up  At 


where 


(3.28) 


« W - — 

and^»  is  the  vertical  displacement. 


For  an  airfoil  undergoing  flapping  and  pitching  oscillations  about 


mid-chord  position 


/ JiioT  , JLijT 

= 26  + <»<  X e 


(3.29) 


Equations  (3.28)  and  (3.29)  are  combined  to  yield 

s ^ + (^1  + J — 


(3.30) 


where  Z'  and  a'  are  the  amplitudes  in  flapping  and  pitching  motions. 


respectively. 

For  oscillatory  (gust)  flow 


6 = x'  e 


(3.31) 


3 

1 

1 
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Equations  (3.28)  and  (3.31)  are  combined  to  yield 


U/=  p<  e e 


//« 


(3.32) 


Lift  and  Moment  Coefficients 


and 


k = = U£Ke 

^ + u i 2’c*; 


i (>A + wt;) 


where  a(x)  is  the  pressure  difference. 

Equations  (3.33)  and  (3.34)  are  combined  to  yield 

- . .2. 


(3.33) 


(3.34) 


Where 


•iAX 


(3.35) 


K=  ice 

Then  the  lift  and  moment  coefficients  ere  evaluated  using  equations 
(3.15)  and  (3.16),  respectively.  Some  results  and  computational  times 
are  presented  below  in  tabular  form  for  both  cases  of  oscillating  airfoils 

* t 

and  oscillatory  flows. 
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Oscillating  Airfoils 


(li 

n 

^mz 

^ma 

0.1 

0.01 

0. 073+0. 525i 

5.297-0.4911 

0.058+0.2861 

2.888-0. 602i 

0.1 

0.50 

0.116+0.5691 

5.760-0.9011 

0.088+0.3071 

3.114-0.9231 

0.1 

0.60 

0. 141+0. 591i 

5.996-1.1441 

0.107+0.3161 

3.219-1.1211 

0.1 

0.70 

0.179+0.6181 

6.296-1.5191 

0.136+0.3261 

3.332-1.4331 

0.1 

0.80 

0.242+0.6501 

6.652-2.1441 

0.184+0.3301 

3.392-1.9671 

0.2 

0.01 

0.102+0.9241 

4.730-0.0881 

0.125+0.5011 

2.582-0.7091 

0.2 

0.50 

0.183+0.9771 

5.050-0.4831 

0.193+0.5171 

2.695-1.0991 

0.2 

0.60 

0.225+0.9981 

5.189-0.6901 

0.231+0.5191 

2.718-1.3231 

0.2 

0.70 

0.288+1.0251 

5.372-1.0081 

0.288+0.5121 

2.706-1.6691 

0.2 

0.80 

0.384+1.0461 

5.550-1.5221 

0.373+0.4731 

2.519-2.2241 

0.3 

0.01 

0.040+1.2771 

4.416+0.4471 

0.176+0.6901 

2.422-0.7411 

0.3 

0.50 

0.154+1.3401 

4.721+0.0731 

0.286+0.6951 

2.488-1.1971 

0.3 

0.60 

0.214+1.3711 

4.879-0.1351 

0.347+0.6871 

2.485-1.4641 

0.3 

0.70 

0.302+1.4031 

5.061-0.4541 

0.435+0.6551 

2.396-1.8691 

0.3 

0.80 

0.444+1.4181 

5.218-1.0201 

0.554+0. 548i 

1.996-2.4841 

0.4 

0.01 

-0.112+1.6011 

4.201+1.009i 

0.210+0.8621 

2.325-0.7571 

0.4 

0.50 

0.041+1.7041 

4.602+0.6261 

0.374+0.8621 

2.394-1.2991 

0.4 

0.60 

0.125+1.7511 

4.802+0.3991, 

0.465+0.8411 

2.366-1.6231 

0.4 

0.70 

0.254+1.7981 

5.032+0.0161 

0.592+0.7731 

2.196-2.1121 

0.4 

0.80 

0.468+1.7941 

5.149-0.7001 

0.726+0.5631 

1.522-2.7471 

0.5 

0.01 

-0.345+1.9241 

4.081+1.5641 

0.233+1.0331 

2.288-0.7751 

0.5 

0.50 

-0.142+2.0891 

4.619+1.1551 

0.465+1.0281 

2.367-1.4231 

0.5 

0.60 

-0.025+2.1601 

4.876+0.8831 

0.593+0.9861 

2.306-1.8171 

0.5 

0.70 

0.171+2.2261 

5.164+0.3831 

0,765+0.8611 

2.021-2.3991 

0.5 

0.80 

0.477+2.1701 

5.160-0.5351 

0.870+0.5231 

1.030-2.9371 

The  computational  time  for  all  the  above  25  cases  was  4.39  seconds 
using  FORTRAN  G Compiler.  In  this  case,  the  airfoil  is  divided 
into  11  strips. 


Oscinatory  Flow 


M 

0.10 

0.01 

0.10 

0.50 

0.10 

0.60 

0.10 

0.70 

0.10 

0.80 

0.50 

0.01 

0.50 

0.50 

0.50 

0.60 

0.50 

0.70 

0.50 

0.80 

1.00 

0.01 

1.00 

0.50 

1.00 

0.60 

1.00 

0.70 

1.00 

0.80 

5.188-0.9671 
5.597-1.4101 
5.800-1.6681 
6.049-2.0571 
6.319-2.6821 
3.391-0.1591 
3.520-0.5261 
3.554-0.7111 
3.545-0.9731 
3.367-1 .2801 
2.404+0.9241 
2.736+0.4591 
2.742+0.1831 
2.568-0.0641 
2.462-0.1101 


2.812-0.5491 

3.007-0.8381 

3.093-1.0101 

3.177-1.2771 

3.196-1.7111 

1.794-0.1421 

1.716-0.4081 

1.629-0.5261 

1.441-0.6431 

1.095-0.5761 

1.226+0.4511 

1.126+0.1661 

0.974+0.1021 

0.784+0.2001 

0.841+0.3541 


The  computational  time  for  all  of  the  above  15  cases,  taking  11 

• i 

strips  over  the  airfoil,  was  4.09  seconds. 

The  present  lifting  surface  method  resulted  In  the  development  of 
a very  efficient  computer  program  for  a 2-d1mens1onal  compressible  flow 
case,  which  In  general  has  one  of  the  most  complicated  formulations. 


IV.  ROTOR  APPLICATIONS 

i)  Steady  State  Application  to  Marine  Propellers 

Elchuri®** applied  the  velocity  potential  formulation  and  the  as- 
sociated numerical  technique  to  the  problem  of  estimating  the  steady  and 
unsteady  loading  distribution  on  typical  ship  propellers.  In  order  to 
study  the  response  of  the  method,  two  B series  propellers,  B-3.35  and 
B-3.50,  were  examined.  Expanded  outlines  of  these  two  types  of  blades 
are  shown  in  Figs.  2 and  3 . As  a first  approximation,  blades  were 
treated  as  thin  lifting  surfaces  with  no  wake.  The  lifting  surface  was 
divided  into  a grid  of  100  boxes  with  ten  such  boxes  along  each  chord 
at  an  equal  number  of  s panwise  locations.  To  obtain  the  distribution 
of  the  discontinuity  in  tne  velocity  potential  on  the  blade  surface, 
the  presence  of  all  the  blades  and  their  wakes  was  considered.  The  doublet 
distribution  obtained  as  a result  of  solving  a set  of  simultaneous 
equations  as  d scussed  earlier  is  related  to  the  differential  pressure 
distribution  on  the  blade  at  any  point  P as 

•f  — / «>  ' \ > 

where  the  relative  free  stream  velocity  is 

^o.  = ^ + . 

The  thurst  produced  by  a blade  section  at  a mean  radius  r^^^  is  the  component 
of  the  resultant  of  lift  and  drag  forces  along  the  direction  of  propeller 
advance.  Therefore,  per  unit  span, 

T = s 

const.  ^ » 

where  the  lift  per  unit  span  due  to  a box  with  P as  the  collocation  point  is 

I ^ [w„^ki-  w.,k,  -ik  ^ ds] 
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and  the  profile  drag  per  unit  span  at  P is 

with  the  local  drag  coefficient  generally  stated  as  an  empirical  function 
of  the  angle  of  attack. 

Torque  necessary  to  overcome  water  resistance  per  unit  span  to  turn 
the  propeller  at  a constant  angular  velocity  is 

5=2:  [•t-VoxCL+-B;}»k]k 

COMSt. 

where  is  the  position  vector  of  point  P with  reference  to  origin  0. 

It  is  customary  to  describe  the  propeller  characteristics  in  terms 
of  thrust  and  torque  coefficients  (aggregated  over  the  entire  span),  and 
the  propeller  efficiency,  as  functions  of  the  advance  coefficient  J. 

These  parameters  are  defined  as. 


Thrust  coefficient,  “ 

Torque  coefficient,  r 


Propeller  efficiency. 


2:  a/ 

TT.  Va  _ Kr  J 


where  the  advance  coefficient  is 

J=  Y»/Cnp)  . 

Thrust  and  orque  coefficients  together  with  propeller  efficiency  as 
functions  of  the  advance  coefficient  at  a pitch  to  diameter  ratio  of  1.0 
for  B-3.35  and  B-3.50  marine  propellers  are  shown  in  Figs.  4 and  5. 

A good  comparison  with  the  experimental  results  of  reference  54  reflect 
upon  the  feasibility  of  the  present  lifting  surface  method  with  regard 
to  its  applicability  to  marine  propellers. 
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ii)  Analysis  of  Unsteady  Airloads  of  Helicopter  Rotors  in  Hover 

Schatzle®®  applied  the  present  numerical  lifting  method  to 
predict  the  unsteady  airloads  on  a multi-bladed  helicopter  rotor 
in  hover  using  a realistic  wake  representation.  The  method 
employed  a velocity  potential  formulation  and  the  wake  geometry 
was  prescribed  using  Landgrebe's  model.  A summary  of  the  results 
were  presented  at  the  AIAA  15th  Aerospace  Sciences  Meeting  (Ref.  56, 
also  to  be  published  in  the  Journal  of  Aircraft) . A brief  discussion 
of  the  method  and  the  results  are  presented  in  this  section. 

Basic  Equations  and  Numerical  Method 

The  relation  between  the  downwash  velocity  at  a point  on  a 
rotor  blade  (wing)  and  distribution  of  doublets  over  the  blade  and 
wake  surfaces  is  derived  in  a previous  section  (Eq.  2.23)  and  is 
given  by 


•h 

h/a4U. 


where  K,  the  local  doublet  intensity,  is  equal  to  the  discontinuity 
in  the  transformed  potential  across  the  wing  or  wake.  The  value 
of  K at  a point  in  the  wake  may  be  determined  from  the  value  of 
K at  the  trailing  edge  of  the  wing  in  the  following  manner.  From 
Euler's  equation,  the  local  pressure  difference  across  a thin 
surface  is  related  to  the  corresponding  discontinuity  in 
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potential  by 

Since  the  wake  cannot  support  a pressure  discontinuity,  =0  in 
eq.  (4.2)  and  thus. 


K(x,y)  = 


In  the  numerical  technique  developed,  it  is  assumed  that 
there  are  M panels  along  the  chord  of  the  blade  and  N panels 
along  the  span.  For  a panel  mn  the  local  doublet  intensity, 
is  constant  over  the  panel  and  may  therefore  be  removed 
from  the  integration  of  eq.  (4.4)  over  that  panel.  Similarly, 
Kte  ^ be  expressed  as  and  removed  from  the  integration 

of  eq.  (4.4)  over  that  strip.  If  the  integrations  in  eq.  (4.4) 
are  carried  out  over  each  panel  and  strip  and  then  combined  to 


form  the  entire  surface  integration,  eq.  (4.4)  may  be  written 
in  influence  coefficient  form  as 
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N 


N 


4irw^-  5, V*'-.  “ “ 


where 


2>Z 


\ J (4.5a) 


and 


(4.5b) 


Here  represents  the  downwash  velocity  induced  at  P due  to 

a unit  velocity  doublet  at  the  panel  mn,  while  ^ represents 
the  downwash  velocity  at  P due  to  a unit  intensity  doublet  in 
the  wake  strip  n.  For  an  M x N panel  lifting  surface,  P takes 
M X N different  values  and  the  problem  reduces  to  a set  of 
simultaneous,  linear  algebraic  equations 


4-Tr\.w]  = 


o 


(4.6) 
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where 


C 


t[^^nnfN*v\ 


^t5  + S s 


nrv^  :;^/V^ 
fYY\  =.  K\ 


The  flow  tangency  condition,  as  detailed  later,  requires  the 
downwash  velocity  induced  by  the  doublets  to  equal  the  downward 
velocity  of  the  wing  at  every  point.  Consequently,  for  prescribed 
wing  motion  and  geometry,  eq.  (4.6)  may  be  solved  for  the  doublet 
distribution,  K. 

The  proper  approach  to  the  solution  of  eq.  (4.6)  is  to  deter- 
mine not  only  the  doublet  distribution  but  the  wake  geometry  as 
well  since  it  is  known  a priori.  This  additional  requirement  can 
increase  computational  effort  appreciably.  However,  for  moderate 
to  high  aspect  ratio  fixed  wings  it  has  been  found  sufficiently 
accurate  to  specify  the  wake  to  be  a planar  vortex  sheet  rather 
than  to  calculate  its  actual  shape.  A similar  approach  is  therefore 
chosen  for  the  rotor  wake.  That  is,  in  the  interest  of  reducing 
computer  time  without  compromising  accuracy,  the  rotor  wake  will 
be  prescribed  rather  than  calculated.  The  difficulty,  however, 
lies  in  choosing  the  proper  wake  shape  for  a given  flight  con- 
dition. The  strip  theory  model  is  simple  but  not  physically  . 
realistic  since  it  ignores  wake  curvature.  The  classical  wake 
model  is  an  improvement,  but  it  neglects  tip  vortex  effects  and 
the  contraction  of  the  wake  of  the  wake  beneath  the  rotor  disc. 

Only  Landgrebe’s  experimentally  determined  rotor  wake  model 
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includes  all  of  these  important  effects.  He  models  the  wake  as 
having  a strong  tip  vortex  and  a weaker  inboard  vortex  sheet 
and  gives  empirical  equations  for  the  radial  and  axial  coordinates 
of  both  as  functions  of  rotor  geometry  and  operating  conditions. 

The  details  of  the  equations  were  presented  in  Ref.  55.  In  this 
report  these  equations  were  used  with  some  slight  modifications. 

The  first  modification  was  to  extend  the  outer  vortex  sheet 
boundary  to  the  tip  vortex,  since  there  is  likely  a physical 
connection  between  the  two.  The  second  was  to  simulate  the  roll 
up  of  the  the  vortex  sheet  into  a tip  vortex  by  gradually  combining 
the  outer  vortex  sheet  wake  strips  into  a single  strip  with  the 
combined  strength  of  the  composite  strips.  This  process  was  done 
over  an  arbitrarily  chosen  one-sixth  rotor  revolution  to  insure 
that  the  tip  vortex  reached  full  strength  before  passing  beneath 
the  following  blade.  The  portion  of  the  wake  which  combined  to 
form  the  tip  vortex  was  determined  by  the  spanwise  location  of 
the  peak  load  on  the  rotor  blade,  and  since  the  peak  load  is  gene- 
rally near  the  85-90%  radius  position,  all  wake  strips  outboard 
of  that  location  were  blended  into  the  tip  vortex.  Furthermore, 
the  outermost  segment  of  the  strips  in  this  region  was  assigned 
tip  vortex  coordinates,  while  all  others  were  given  vortex  sheet 
coordinates.  It  was  also  assumed  that  the  vortex  sheet  contracted 
radially  at  the  same  rate  as  the  tip  vortex. 

The  blade  surface  is  divided  into  a number  of  small  panels 
and  the  panel  size  is  smaller  in  regions  where  K varies  most 
rapidly,  i.e.,  near  the  blade  edges.  The  freestream  velocity  is 
assumed  to  be  constant  across  a panel  and  equal  to  the  velocity 
at  the  geometric  centroid  of  the  panel.  This  velocity  varies  with 
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span  due  to  the  nature  of  the  rotor  flow  field.  Extending  from 
each  trailing  edge  panel  is  a wake  strip  whose  coordinates  are 
determined  by  Landgrebe's  equations  (Fig.  6 shows  the  panel 
approximation  of  the  rotor  blade  and  the  modified  wake  model) . 

Calculation  of  Influence  Coefficients  j 

1 

Eq.  (4.6)  can  be  solved  for  K at  a discrete  number  of  points  \ 

on  the  surface  which,  placed  at  the  geometric  centroid  of  each  panel. 

The  influence  of  a surface  panel  on  a collocation  point  is  given 
by  eq.  (4.5a),  which  is  evaluated  numerically  using  a two-dimensional 
Gaussian  quadrature.  The  influence  of  a wake  strip  on  a collocation 
point  is  given  by  eq.  (4.5b)  with  the  following  modification:  implicit 
in  the  use  of  (4.3)  for  enforcing  the  zero-pressure  condition  in 
the  wake  is  the  assumption  that  the  local  streamlines  are  in  the 
direction  of  the  positive  X-axis.  In  a linearized  sense  this  is 
true  for  fixed  wings.  However,  for  rotary  wings  the  streamlines 
are  directed  roughly  helical  paths  beneath  the  rotor  disc.  To 
reflect  this  difference  eq.  (4.3)  is  therefore  modified  to 

kCsi  = I 

where  S is  the  distance  measured  along  the  centerline  of  any  | 

wake  strip.  Of  course,  eq.  (4.5b)  must  be  changed  in  accordance  : 

with  eq.(4.7).  The  effect  that  a particular  wake  strip  has  on 
a collocation  point  is  therefore  computed  by  adding  the  effect 
that  each  panel  in  the  strip  has  on  that  point  using  the  modi- 
fied form  of  eq.(4.5b).  As  with  the  surface  panels,  the  inte- 
grations are  carried  out  with  a two-dimensional  Gaussian  quad- 


rature . 
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Surface  Boundary  Condition 

The  vertical  displacement  of  the  blade  is  assumed  to  have 
a steady-state  component  due  to  coning  and  angle  of  attack,  and 
an  unsteady  component  due  to  flapping  and  twisting  about  this 
stady  position,  i.e.. 


'T  -r  ^ 


(4.8) 


It  can  be  seen  from  Fig.  7 that 


(4.9a) 

- yC^) 

(4.9b) 

tc 

(4.9c) 

The  relations  between  the  downwash  velocity 

Wp  and  "y , and 

and  the  transformed 

downwash  velocity  W are, 

P 

respectively , 

H - 

H-  U ^ 

(4.10) 

“•  L ^ ^ 

= 

Wj.  e 

(4.11) 

Eqs.(4.8)  through  (4.11)  are  combined  to  give 


where 
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The  solution  for  steady  flow  is  obtained  by  settingl=0  andX-0 
in  eqs.(4.5a)  and  (4.5b)  and  using  eqs . (4,6)  and  (4.12a).  Simi- 
larly, the  solutions  for  flapping  and  twisting  are  obtained  by 
using  eqs. (4.12a)  and  (4.12b),  respectively  in  eq.(4.6). 


Aerodynamic  Loads 

Eq.(4.2),  the  linearized  expression  for  the  lift  in  trans- 


formed coordinates,  is  modified  to 


^Cx.YjT) 


\ i (ax  -i-uST ) 


(4.13) 


By  integrating  eq.(4.13)  along  the  chord,  the  local  lift  and 


moment  per  unit  span  may  be  expressed  as  I 


tto  1 


(4.14) 
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The  unsteady  lift  and  moment  coefficients  for  flappinf  and 
twisting  motions  are  expressed  as 

iJ  ' 

pjjipi  ~ ^nz^fctr  , (4.17) 

where  z'  . =2'  . //.  It  should  be  noted  that  C_  , C_  , C,,  , 

tip  tip  t Lz  Mz 

and  are  complex  quantities . 

Solution  Procedure 

For  a given  rotor  geometry  and  flight  condition,  blade 
element  theory  is  used  to  calculate  an  initial  guess  to  the  rotor 
thrust  coefficients,  C^.  Based  on  this  value  of  and  the 
initial  conditions  an  initial  wake  shape  is  prescribed  using 
the  modified  Landgrebe  model.  The  resultant  steady-state  airloads 
are  then  calculated  and  an  improved  guess  for  C,j,  is  made.  If  this 
second  value  differs  appreciably  from  the  first,  a new  wake 
geometry  exists  which  is  compatible  with  the  steady-state  load 
distribution,  i.e.  , The  unsteady  aerodynamic  derivatives 

are  then  calculated  based  on  this  converged  geometry. 

Results 

Steady  airloads  and  unsteady  aerodynamic  derivatives  were 
calculated  for  an  XH-51A  helicopter  rotor  in  hover.  This  rotor 
was  chosen  because  experimental  data®^  for  the  steady-state 
load  distribution  was  available  for  comparison  with  theory. 
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The  blade  geometry  and  flight  conditions  were: 


NB 

Pi 

II 

= 17.5’ 

c 

= 1.08' , y 

' CO 

= 2.33' 

01 

CD 

0 

in 

II 

= 10.61' 

Nr 

= 355  (M^,p  = 

0.58)  . 

In  the  unsteady  case  the  blades  were  assxamed  to  undergo  rigid 
flapping  and  torsional  motion  such  that 

fO  , y -s  e 

(y-e)  , y e 

(R-e) 

and 

F(y)  = 1. 

Figs.  8 and  9 show  the  calculated  lift  and  moment  distributions 
for  the  XH-51A  using  the  present  theory  with  the  classical  and 
Landgrebe  wake  models.  Both  models  overpredict  the  overall  blade 
lift  although  Landgrebe 's  model  predicts  the  loading  trends 
more  accurately,  specifically  the  load  spike  near  the  tip  due  to 
the  tip  vortex.  The  mean  value  for  for  the  classical  model 
was  0.00577  and  that  for  the  Landgrebe  model  was  0.00584,  while 
the  experimental  value  was  roughly  0.005  - - a difference  of  15 
to  17%  between  theory  and  experiment.  This  difference  may  be  due 
to  the  approximate  nature  of  the  present  method,  i.e.,  discretizing 
the  linear  spanwise  velocity  field  into  constant  velocity 
segments.  No  comparison  was  made  between  the  theoretical  and 
experimental  moment  distributions  since  the  measured  data  were 
referenced  to  an  unspecified  axis.  It  might  be  expected. 
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Figure  9.  Moment  per  Unit  Span,  XH-51A  Helicopter  Rotor 
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however,  that  Landgrebe's  model  would  give  better  results  in  the 
tip  region  than  would  be  the  classical  model.  The  number  of  wake 
revolutions  used  in  this  study  was  varied  from  one  to  four  with 
no  appreciable  change  in  the  overall  thrust  or  character  of  the 
loading  curves,  but  with  a drastic  increase  in  computational  time. 
The  displyed  results  are  therefore  those  obtained  using  a single 
wake  revolution.  For  lightly  loaded  rotors  (where  the  rotor  wake 
is  closer  to  the  blade)  it  may  be  necessary  to  include  more  wake 
revolutions  to  get  accurate  results.  The  calculation  of  the  steady 
state  loads  using  either  wake  model  required  approximately  36 
seconds  of  CPU  time  in  G-Level  FORTRAN  on  Amdahl  470V/6  computer. 

Based  on  the  steady-state  results  only  the  Landgrebe  model 
was  used  for  three  frequency  ratios  shown  in  Figs. 10  through  17. 
Evident  in  each  of  these  figures  is  the  strong  tip  vortex  influence 
on  the  loading  over  the  outer  20%  of  the  blade.  In  the  present 
study  no  attempt  is  made  to  explain  the  variation  of  these  quanti- 
ties as  a function  of  frequency  ratio  or  its  significance  in  flutter 
analysis.  However,  based  on  the  reasonable  correlation  between 
steady-state  theoretical  and  experimental  loading  using  the  Land- 
grebe  model,  it  is  hoped  that  the  computed  unsteady  aerodynamic 
derivatives  will  also  be  reasonably  accurate.  In  the  future 
these  results  will  be  compared  with  those  of  strip  theory  and 
classical  wake  model  representations  and  applied  to  investigate 
their  effect  on  aeroelastic  stability  boundaries  as  was  done  in 
Ref.  58.  The  calculation  of  the  steady-state  airloads  and  wake 
geometry  and  the  resulting  unsteady  aerodynamic  derivatives 
required  approximately  2.8  minutes  of  CPU  time  in  G-Level  FORTRAN 


Figure  15.  Imaginary  Part  of  Moment  Due  to  Flapping 
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on  an  Amdahl  470V/6  computer. 

Conclusions 

A velocity  potential  lifting  surface  method  has  been  used 
in  conjunction  with  a realistic  rotor  wake  model  to  compute  the 
steady  airloads  and  unsteady  aerodynamic  derivatives  for  an 
arbitrary  hovering  rotor  in  compressible  subsonic  flow.  The 
theoretical  results  compare  well  with  experimental  data  for  four- 
bladed  rotors,  suggesting  that  unsteady  loads  calculated  with  the 
same  wake  geometry  may  also  be  reasonably  accurate.  With  the 
proper  choice  of  flappinf  and  torsional  modes,  the  present  method 
could  be  used  to  generate  aerodynamic  derivatives  for  use  in 
flutter  analysis.  In  fact  the  next  section  deals  with  this  topic. 

iii)  Aeroelastic  Analysis  of  a Sinqle-Bladed  Rotor 

The  unsteady  aerodynamic  analysis  and  the  computer  program 
described  in  the  previous  section  is  applied  to  predict  the  flutter 
speed  of  a single-blades  rotor  in  Ref.  59.  The  results  an!  the 
discussion  of  this  analysis  is  presented  in  this  section. 

A single-bladed  rotor  with  the  following  geometric  and  stru- 
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(radius  of  gyration  ) 
(bending/ torsion  freq.  ratio) 
(density  ratio) 

(aft  of  elastic  axis) 


ctural  properties  was  considered: 
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To  assure  compatibility  between  the  two  and  three-dimensional 
aerodynamic  derivatives,  a classical  wake  was  prescribed  in  the 
three  dimensional  analysis.  The  use  of  classical  wake  resulted 
in  no  loss  in  generality  of  the  three-dimensional  theory  since 
it  has  been  shown  that  there  is  little  difference  between  the 
results  obtained  from  using  either  Landgrebe's  wake  model  or 
a classical  wake  (Ref.  59)  . 

The  principal  flutter  parameters  in  this  analysis  were  the 
location  of  the  elastic  axis  (a)  and  the  bending/torsion  frequency 
ratio  (^) . The  varation  of  flutter  speed  is  shown  in  Fig.  18  for 
the  assumed  blade  section  at  0.75  R,  using  both  two-  and  three- 
dimensional  aerodynamics.  Using  two-dimensional  aerodynamics 
for  locations  of  the  elastic  axis  ahead  of  the  quater  chord,  an 
analyst  could  mistakenly  conclude  that  there  is  no  chance  for 
flutter  to  occur;  for  this  case  the  three-dimensional  aerodynamics 
indicates  a lower  flutter  speed.  However,  when  the  elastic 
axis  is  moved  aft  of  a = 0.3,  three-dimensional  aerodynamics  is 
less  conservative  than  two-dimensional  aerodynamics.  Consequently, 
two-dimensional  theory  should  be  used. 

For  the  rest  of  this  study,  the  elastic  axis  was  located  at 
a = 0.2.  Fig.  19  shows  the  variation  of  the  flutter  speed  with 
the  ratio  of  the  bending/torsion  frequencies.  It  is  seen  that 
two-dimensional  theory  predicts  a lower  flutter  speed  than  does 
the  three-dimensional  theory  for  all  values  of  C.  Also,  for 
values  of  ^ near  1.0  and  greater,  both  aerodynamic  theories 
predict  no  flutter  for  values  of  reduced  frequecy  less  than  0.5 
These  same  trends  are  also  seen  when  the  flutter  frequency  is 
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plotted  against  the  bending/torsion  frequency  ratio  (Fig.  20) . 

In  terms  of  a rotor  system,  these  results  indicate  that  no 
flutter  condition  will  be  encountered  for  the  bending  stiffness 
greater  than  the  torsional  rigidity.  Further  from  these  two 
figures,  one  can  conclude  that  the  finite  aspect  ratio  effects 
are  such  that  they  reduce  the  flutter  speeds. 

Fig.  21  shows  the  variation  of  the  flutter  speed  with  the 
density  ratio.  For  both  cases,  the  flutter  speed  increases  with 
an  increase  in  the  density  ratio,  and  the  gradient  for  the  two- 
dimensional  curve  is  seen  to  be  decreasing  faster  than  that  of  the 
three-dimensional  curve.  The  trend  indicates  that  the  increase 
in  density  ratio  has  a slightly  less  stabilizing  effect  on  the 
airfoil  motion  when  using  two-dimensional  airfoil  theory.  Thus, 
for  high  values  of  density  ratio,  a very  conservative  flutter 
speed  may  be  predicted  using  strip  theory. 

The  varaition  of  the  flutter  speed  with  the  square  of  the 
non-dimensional  radius  of  gyration  is  shown  in  Fig.  22.  Once 
again,  the  flutter  speed  obtained  using  strip  theory  is  the 
conservative  result.  The  flutter  speed  increses  as  the  radius 
of  gyration  increases. 

The  varaition  of  c.g.  location  on  the  flutter  speed  is 
shown  in  Fig.  23.  Here,  using  three-dimensional  aerodynamic 
derivatives,  there  is  insignificant  variation  in  the  flutter 
speed.  However,  for  strip  theory,  the  effect  of  c.g.  location 
variation  has  a much  more  pronounced  effect  on  the  flutter 
speed. 

The  results  presented  in  this  section  are  very  preliminary 
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